Read in data

Subsets only the NK and T cells.

library(Seurat)
## Attaching SeuratObject
## Attaching sp
library(dittoSeq)
## Loading required package: ggplot2
library(readxl)
source("../functions_dge.R") |> suppressPackageStartupMessages()
## Warning: package 'scuttle' was built under R version 4.2.1
## Warning: package 'IRanges' was built under R version 4.2.1
## Warning: package 'GenomeInfoDb' was built under R version 4.2.1
## Warning: package 'scran' was built under R version 4.2.1
## Warning: package 'edgeR' was built under R version 4.2.1
## Warning: package 'limma' was built under R version 4.2.1
sc <- readRDS("../analysis/sc_integrated_annotated.rds")
tc <- subset(sc, subset = celltype_fine %in% c("NK cells", "T cells"))

rm(sc)

# load in file that can be used to re-orde samples in e.g. barplot
sample_order <- read_excel("../metadata/sample_order.xlsx")
sample_order$sample <- sapply(sample_order$`Sample ID`, function(x) {
  strsplit(x, " ")[[1]][1]
})

Re-do dimensiontality reduction on cycling cells only

DefaultAssay(tc) <- "integrated"
tc <- ScaleData(tc, verbose = FALSE)
tc <- RunPCA(tc, verbose = FALSE)
tc <- RunUMAP(tc, dims = 1:50)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

UMAPS

library(dittoSeq)
dittoDimPlot(tc, var = "orig.ident", size = 0.25)

dittoDimPlot(tc, var = "Phase", size = 0.25)

dittoDimPlot(tc, var = "celltype_fine", size = 0.25)

Clustering

tc <- Seurat::FindNeighbors(tc, dims = 1:50, reduction = "pca")
tc <- Seurat::FindClusters(tc, resolution = 0.8)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 35543
## Number of edges: 1811733
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8473
## Number of communities: 26
## Elapsed time: 11 seconds
pdf("../plots/umap_t_cells_cl0.8.pdf", width = 9)
p <- dittoDimPlot(tc,
  var = "integrated_snn_res.0.8",
  size = 0.5, do.label = TRUE
)
print(p)
dev.off()
## quartz_off_screen 
##                 2
print(p)

tc$tcell_annotation <- "other"
DefaultAssay(tc) <- "RNA"

Plot marker genes

markers <- read.csv("../metadata/marker_genes_man_tcell.csv")
mli <- split(markers, markers$tcell_annotation)

for (cellt in names(mli)) {
  marker_genes <- mli[[cellt]]$gene
  cellt_name <- gsub("[ \\+]", "", cellt)

  pdf(paste0("../plots/UMAP_tcells_markers_", cellt_name, ".pdf"))
  multi_dittoDimPlot(tc,
    var = marker_genes,
    min.color = "grey",
    max.color = "purple",
    size = 0.25,
    order = "increasing"
  )
  dev.off()
}
dittoDotPlot(tc,
  vars = unique(markers$gene),
  group.by = "integrated_snn_res.0.8"
)

Assign clusters to T cell type

tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(3, 12)] <- "NK cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% 10] <- "T regulatory cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% 9] <- "T effector cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(5, 8)] <- "Naive T cells"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(1, 7, 15, 4)] <- "Cytotoxic T cells CD8+"
tc$tcell_annotation[tc$integrated_snn_res.0.8 %in% c(2, 13)] <- "Exhausted T cells CD8+"
saveRDS(tc, "../analysis/tcells_annotated.rds")

Visualization annotation

p <- dittoDimPlot(tc,
  var = "tcell_annotation", do.label = TRUE, size = 0.25,
  labels.size = 2
)
print(p)

pdf("../plots/umap_manual_annotation_tcells_labeled.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen 
##                 2
p <- dittoDimPlot(tc, var = "tcell_annotation", size = 0.25)
print(p)

pdf("../plots/umap_manual_annotation_tcells.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen 
##                 2
p <- dittoBarPlot(tc,
  var = "tcell_annotation", group.by = "orig.ident",
  x.reorder = match(
    sample_order$sample,
    unique(tc$orig.ident)
  )
)
print(p)

pdf("../plots/barplot_manual_annotation_tcells.pdf", width = 10)
print(p)
dev.off()
## quartz_off_screen 
##                 2
library(dittoSeq)

p <- dittoDotPlot(tc,
  vars = unique(markers$gene),
  group.by = "tcell_annotation",
  min.color = "lightgrey", max.color = "darkred"
)
print(p)

pdf("../plots/dotplot_tcell_markers.pdf")
print(p)
dev.off()
## quartz_off_screen 
##                 2
tmark <- read.delim("../metadata/tcell_markers.txt", header = F)

cr <- colorRampPalette(c("white", "black", "red"))

p <- dittoHeatmap(tc, unique(c(markers$gene, tmark$V1)),
  scaled.to.max = TRUE,
  heatmap.colors.max.scaled = cr(25),
  annot.by = c("tcell_annotation", "integrated_snn_res.0.8")
)

pdf("../plots/heatmap_tcell_markers.pdf", height = 10, width = 15)
print(p)
dev.off()
## quartz_off_screen 
##                 2